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1.  FOREWORD 


The  project  entitled  ”  Large  Deformation  Analysis  of  Nonlinear  Homogeneous 
and  Heterogeneous  Media  Using  an  Arbitrary  Lagrangian-Eulerian  Finite 
Element  Method  "  (Grant  No.  DAAL03-9 1-0-0069)  has  been  transferred  from  the 
University  of  Alabama  to  the  Ohio  State  University  effective  August  16. 1991.  The  principal 
investigatOi  of  the  project.  Dr.  Somnath  Ghosh  has  assumed  a  faculty  position  at  the  Ohio 
State  University  and  the  project  is  continuing  at  OSU.  Therefore,  in  this  report  only  the  first 
part  of  the  ongoing  work  is  presented.  Subsequent  progresses  will  be  reported  from  the 
Ohio  State  University. 

2.  LIST  OF  APPENDICES 

First  draft  of  a  manuscript  entitled "  ANALYSIS  OF  RANDOM  COMPOSITES  USING 
VORONOI  CELL  FINITE  ELEMENTS” 


3.  BODY  OF  THE  PROBLEM  STUDIED 

A.  STATEMENT  OF  THE  PROBLEM 

In  this  part  of  the  study,  a  new  finite  element  technique  has  been  developed  for  macroscopic 
analysis  of  composites  with  a  stochastic  dispersion  of  microscopic  heterogeneities.  In  these 
composites,  particles  or  fibers  of  various  shapes  or  sizes  may  be  randomly  scattered  within 
the  matrix.  The  present  work  introduces  a  material  based  discretization  system  that 
adequately  accounts  for  the  morphology  of  the  heterogeneous  domain  and  then  formulates  a 
finite  element  scheme  for  accommodating  multi-sided  elements  with  second  phase  materials 
present 


B.  SUMMARY  OF  MOST  IMPORTANT  RESULTS 


The  heterogeneous  domain  with  random  distribution  of  second  phase  has  been  discretized 
by  Dirichlet  tessellation  so  that  each  element  contains  only  one  inclusion.  The  resulting 
networic  of  Voronoi  cells  are  now  assumed  to  be  elements  in  a  finite  element  formulation. 

Indeed,  the  multi-noded  Voronoi  cells  are  rather  unconventional  elements  and  a 
displacement  based  formulation  with  such  elements  is  extremely  difficult  due  to  the  varying 
number  of  sides  in  each  element  An  assumed  stress  hybrid  finite  element  formulation  has 
proven  to  be  extremely  successful  in  this  process.  Several  different  test  cases  have  been 
experimented  with  for  linear  elastic  problems  and  good  correlation  between  numerical 
solutions  and  analytical  results  have  been  obtained. 

The  effect  of  the  second  phase  has  been  incorporated  by  the  introduction  of  a 
transformation  strain  in  the  element  level  in  a  fashion  similar  to  self-consistent  schemes. 
The  results  have  been  compared  with  conventional  finite  element  predictions  and  are  in 
reasonable  agreement  Work  is  now  in  progress  for  refining  this  model. 


C.  PUBLICATIONS 

Two  manuscripts  arc  presently  being  written  and  will  be  forwarded  to  ARO  as  soon  as  they 
are  completed. 


D.  DEGREES  AWARDED  FROM  THIS  PROJECT 
No  degrees  have  been  awarded  so  far  from  this  project 

4.  REPORT  OF  INVENTIONS 


A  New  Finite  clement  formulation  for  random  composites  using  Voronoi  Cell  elements 


□  □ 
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1  Introduction 


Research  related  to  the  study  of  composite  behavior  during  the  last  two  decades,  has  r^ulted 
in  a  number  of  computational  models  for  their  analysis.  Most  of  these  models  use  the  homog¬ 
enization  method  for  reflecting  the  influence  of  materi2d  microstructure  on  the  macroscopic 
behavior.  These  methods  are  based  on  asymptotic  analyses,  which  predominantly  make  as¬ 
sumptions  of  periodic  repetition  of  microstructures.  However,  in  practice,  there  exists  a  wide 
variety  of  composite  szmiples  with  a  stochastic  dispersion  of  microscopic  heterogeneities.  In 
thest  composites,  particles  or  fibers  of  \-arious  shapes  and  sizes  may  be  randomly  scattered 
within  the  matrix  and  even  be  clustered  in  certain  regions.  The  distribution  of  shapes,  sizes 
and  the  spatial  coordinates  of  the  second  phase  has  a  profoimd  influence  on  the  mechanical 
beha\*ior  of  the  overall  structure  and  therefore,  must  be  taken  into  account  in  a  rigorous 
anah^is. 

A  fundamental  requirement  in  the  development  of  finite  element  models  for  an2dysis  of 
the  above  class  of  random  composites  is  the  generation  of  a  robust  mesh,  that  will  ad¬ 
equately  account  for  the  morphology  of  the  composite  domain  under  consideration.  The 
Dirichlet  tessellation  seems  to  provide  an  excellent  foundation  for  natural  evolution  of  dis¬ 
cretization  process  while  accoimting  for  the  microstructure.  This  is  a  method  of  subdividing 
an  Euclidean  space  into  n-dimensional  boimded  convex  polytopes.  It  may  be  perceived  as 
production  of  a  network  of  interfaces,  formed  by  the  impingement  of  expanding  hyperspheres 
about  random  nuclei  that  are  growing  at  a  uniform  rate  from  zero.  If  the  second  phase  parti¬ 
cles  are  realized  as  points  in  space,  the  convex  polytopes  (polygons  in  two  dimensions)  known 
as  Voronoi  polytopes  resulting  from  this  discretization,  would  encompass  one  inclusion  each 
at  most.  These  Voronoi  polygons, in  two  dimension,  are  used  as  elements  in  a  finite  element 
analysis  of  random  composites. 
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2  Formulation  of  Voronoi  Cell  Finite  Elements  by 
Medlett’s  Method 

Mallett  [1]  proposed  a  least  square  finite  element  scheme  and  a  finite  difference  element 
scheme  for  the  formulation  of  voronoi  polygon  elements.  In  the  first  scheme,  the  displacement 
components  in  each  element  is  represented  as  linear  function  of  position  and  then  they  are 
determined  by  giving  a  least  square  fit  to  the  nodal  values  of  the  displacement  components. 
The  latter  scheme  completely  eliminates  the  need  to  assume  any  sort  of  approximation  for  the 
displacement  field  in  the  interior  of  the  element.  Instead,  the  element’s  average  displacement 
gradient  components  are  C2Jculated  to  form  the  strain  displacement  matrix. 


2.1  The  Least  Square  Finite  Element  Scheme: 

Displacement  field  in  a  n-sided  polygonal  element  that  represents  the  rigid  body  displace¬ 
ments  and  the  constant  strain  field,  can  be  expressed  as  follows 


+  ^2^  +  oizy 

v{x,y)  -  ^  02X  03y  (1) 

where  x  and  y  are  global  Cartesian  coordinates. 


Evaluating  the  element  displacement  components  at  the  i-th  node  and  reqmring  that 
these  values  coincide  with  the  nodal  displacement  components  u,  and  v,-  leads  to  the  following 
equations 


■ 

1  yi 

.  1  y„ 

Qfi 
^3  . 


or 


{u}  =  [^]{a} 


(2) 


As  the  number  of  nodes  of  an  element  will  be,  in  most,  cases  more  than  three,  equations 
(2)  constitutes  an  overdetermined  set  of  equations.  Applying  least  square  method  to  solve 
the  overdetermined  system,  one  obtains 

{a}  =  ([AnAl)->[A]qu}  (3) 
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Substituting  this  expression  for  {a}  and  the  similar  expression  for  {^}  into  eq.  (1)  gives 


«(*,»)=[!  *  m 

f(x,y)  =  [i  *  y](MnA]nAf {v}  =  iV{v} 


The  components  of  strain  can  be  related  to  nodal  displacement  components  by 


7ry 


r  dN 


^  0 
“  i 


dN 

■5F 


{"} 


(4) 


W  =  (5]{<i}  (5) 

where  [5]  is  the  strain  displacement  matrix  and  {df}  is  he  element  nodal  displacement  vector. 
It  can  be  easily  seen  that  and  is  constant*oveT  the  element  and  so  is  matrix  [B].  The 
element  sti&ess  matrix,  for  unit  thickness,  thus  can  be  expressed  as 


[Kl  =  J[Bm[B]dv  =  MBmiB]  (6) 

where  is  the  2Lrea  of  the  element  and  [J^]  is  the  elasticity  matrix. 

The  nodal  forces  implied  by  the  above  stifeiess  matrix  are  not  statically  equivalent  to 
the  edge  tractions  implied  by  the  stresses  within  the  element.  To  rectify  this  deficiency  a 
[5]  matrix  of  size  2nx3  is  determined  that  transforms  the  stress  components  into  statically 
equivalent  nodal  force  components. 


{F},  =  [5]{<t}  =  [Sll£][B]{<i}  (7) 

from  which  element  stiffness  matrix  can  be  written  as 

[K],  =  [5][£J[5)  (8) 

which  in  general  is  not  symmetric. 


2.2  The  Finite  Difference  Element  Scheme  : 


In  this  scheme,  the  average  of  displacement  gradients  ||,|^  and  |^,|^  over  an  element  is 
calculated  and  written  as 

du  ^  a^u  du  _  b^u 

dx  Ae  dy  Ae 


dv  a^v 
dx  Ae 


dv  ^  b^y 
dy  “  Ae 
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where 


«i  =  \(y{+i  -  Vi-i) 

ii  =  - ij+i)  »  =  l---n 


with  the  underrtiuding  that  (io,y«)  and  {xn+uVn+i)  are, to  be  interpreted  as  (®n,yn)  and 
respectively 

Assuming  a  state  of  constant  strain  within  the  element,  element’s  strain  component 
averages  are  then  given  by 


or 


_L 

A, 


0^ 

6^ 


{e}  =  [B\id} 


(9) 


The  matrix  [S]  is  constant  over  the  element,  hence  the  stiffness  matrix  can  be  expressed 


as 


[K].  =  A,[Bf[E][B] 


(10) 


2.3  Rank  Deficiency  of  Stiffness  Matrix 

The  element  stiffness  matrix  produced  by  either  of  the  schemes  has  a  rank  of  three,  irrespec¬ 
tive  of  the  number  of  nodes  n  in  an  element,  and  its  rank  deficiency  is  2n-6.  Assemblage 
of  rajik  deficient  stiffness  matrices  produces  a  rank  deficient  global  stiffiiess  matrix.  Mallett 
proposed  two  different  methods  for  control  of  hourglass  or  zero  energy  nodes  produced  by 
a  rank  deficient  stiffness  matrix.  In  the  first  method  an  artificial  stiffness  is  added  at  the 
element  stiffness  matrix  level,  which  will  inhibit  hourglass  deformation  without  seriously 
affecting  non-hourglass  deformation.  In  the  second  method  a  constraint  is  imposed  at  the 
global  stiffness  level,  by  means  of  the  Lagrange  multiplier  technique,  to  eliminate  rather 
than  suppress  hourglass  deformation.  However  elimination  of  hourglass  modes  at  the  ele¬ 
ment  level  leads  to  a  global  stiffness  matrix  that  is  excessively  stiff.  Elimination  of  hourglass 
deformation  at  the  global  level  has  been  found  to  be  the  only  viable  alternative  and  the 
following  numerical  examples  are  based  on  the  same  principle. 

2.4  Numerical  Examples 

The  material  properties  assumed  here  are  used  in  all  subsequent  examples  unless  otherwise 
specified.  Further  the  problems  are  always  solved  under  plain  strain  condition.  Poisson’s 
ratio  is  assumed  as  0.3  and  Young’s  modulus  as  2  KN/sq.mm.  The  value  of  Young’s  modulus 
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is  unimportant  as  it  merely  serves  as  a  scale  factor  for  the  displacement. 

Stretching  Problem: 

A  square  material  body  of  size  40  mm  x  40  mm  is  discretized  into  24  voronoi  polygons 
(Fig.  1)  and  is  subjected  to  a  uniform  tensile  force  of  100  KN.  Analysis  of  this  problem  by 
MaUett’s  method  yields  end  deflection  of  4.55  mm  in  X  direction  and  (t*  of  2.5  kN/sq.mm 
in  individual  elements  which  are  the  exact  theoretical  values.  Also  the  displacement  of  the 
nodes  in  X  direction  versus  their  position  from  the  fixed  end  has  been  found  to  be  linearly 
proportional  (Fig.  2). 


>1 


Fig.  1.  Body  subjected  to  tensile  load. 


Fig.  2.  Nodal  displacement  in  X  direction  vs.  nodal  position. 
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Bending  Problem: 

A  material  body  of  size  40  mm  x  160  mm  is  used  as  a  cantilever  beam  with  one  end 
fixed  and  other  end  subjected  to  pure  moment  of  40  KN-mm.  The  beam  is  discretized  into 
28  voronoi  polygon  elements  (Fig.  3).  Longitudinal  stress  Cx  along  a  transverse  section  is 
calculated  by  the  above  method  and  is  plotted  as  shown  in  Fig.  4.  The  stresses  are  found 
to  be  constant  within  an  element  and  differ  widely  from  theoretical  value. 


Fig.  3.  Beam  subjected  to  pure  end  moment. 


Fig.  4.  Sigma  X  along  a  transverse  section. 
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OeflcctiM  la  Y  dlrcctloo 


The  calculated  tip  deflection  is  1.546  mm,  which  is  only  35%  of  the  theoretical  value  of 
4.368  mm.  Fig,  5  shows  a  plot  of  the  vertical  displacement  of  the  nodes  at  the  top  surface 
against  the  theoretical  values. 


Fig.  5.  Deflection  of  top  surface. 
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Proposed  Formulation  by  Assumed  Stress  Hybrid 
Method 

Large  amount  of  error  in  predicting  the  deflection  and  stress  of  a  cantilever  beam  under  pure 
end  moment  by  Mallett’s  method  was  the  motivation  to  look  for  an  alternative  formulation. 
In  the  earlier  section  it  was  shown  that  a  least  square  fit  of  the  displacement  field  from 
the  nodal  values  produce  a  rank  deficient  sti&ess  matrix  for  elements  having  more  than 
three  nodes.  Even  if  the  rank  defficiency  is  taken  care  of  by  global  hourglass  control,  the 
elements  exhibit  constant  stress  which  are  widely  different  from  theoretical  values.  Pian  [2] 
suggested  an  alternative  derivation  of  the  element  stiffness  matrices  by  assumed  stress  hybrid 
method.  In  the  proposed  method,  instead  of  a  required  continuous  displacement  function 
over  the  element,  it  is  necessary  only  to  write  down  the  boundary  displacements  that  will 
guarantee  a  complete  displacement  compatibility.  The  assumed-stress  hybrid  method  is  a 
way  of  formulating  a  sti&ess  matrix  by  use  of  independent  assumptions  of  (a)  an  equilibrium 
stress  field  within  the  element,  and  (b)  interelement-compatible  displacement  nodes  on  the 
element  boundary.  As  this  method  necessitates  the  interpolation  of  displacement  field  not  in 
the  interior  of  the  element  but  only  along  the  element  boundary,  it  becomes  easy  to  handle 
polygonal  elements  having  any  number  of  nodes. 

3.1  Element  Formulation 

The  derivation  is  based  on  the  principle  of  minimum  complementary  energy,  for  which  the 
functional  to  be  varied  is 


Tg  —  /  f  TiUids 

Jv  It  Jdv 


(11) 


where  Sijki  is  the  elastic  compliance,  tJ,*  are  the  prescribed  displacements  along  the  boundary 
dv  and  Ti  are  the  boundary  tractions  that  are  related  to  <Tij. 

The  stresses  satisfy  the  equilibrium  condition 


dcTjj 

dxj 


=  0 


(12) 


and  is  compatible  with  the  prescribed  boundarj'  tractions.  In  applying  the  finite  element 
method,  the  assumed  stress  field  need  not  be  continuous  across  the  interelement  boundaries, 
but  equilibrium  must  be  maintained  for  the  surface  tractions  T,-,  defined  by 


Ti  =  (13) 

where  rij  are  the  components  of  the  unit  vector  normal  to  the  boundary. 

In  the  application  of  this  variational  principle,  one  begins  with  a  stress  field  that  satisfies 
the  differential  equations  of  equilibrium.  The  stress  distribution  can  be  expressed  as 


(14) 
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where  {fi}  is  the  colunm  of  m  imdetenmned  stress  coefficients  h"'  Pm  and  terms  of  the 
matrix  [P]  are  functions  of  coordinates  Xi^yi, 

The  surface  tractions  T,*  have  been  expressed  in  terms  of  the  stress  by  equation  (13) 
and  hence  can  be  related  to  the  undetermined  stress  coefficients  as 

{T}  =  [Rm  (15) 

Also,  the  prescribed  boundary  displacements  Ui  can  be  interpolated  from  generalized 
displacement  {g}  at  the  nodes,  in  the  form 

{«}  =  [i]{9}  (16) 

where  the  terms  in  the  matrix  [Z]  contain  coordinates  on  the  surface. 

The  equation  (11)  can  be  rewritten  as 

= L  (17) 

Substituting  the  expressions  for  <t,T  and  u  one  obtains 

’Tc = -  {mG]U}  (18) 

where  [H]  =  l[Pf[S][P]dv  (19) 

and  [G]  =  (20) 

Making  tTc  stationary  with  respect  to  stress  parameters 

dv 

^  =  0  for  i=l,2,---m  (21) 

gives 

wm  =  [G1{9}  or  {0}  =  [FI-nCHg}  (22) 

Substituting  {0}  in  the  expression  for  complementary  strain  energy  for  an  element  of 
volume  u,  one  obtains 


=  \{0Vm0} 
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where  [K]  =  [Cf^[HY'^[G\  is  the  stifeess  matrix. 

The  stiffness  matrix  [K]  will  be  rank  deficient  if  it’s  rank  is  less  than  n-/  where  n  is  the 
number  of  degrees  of  freedom  and  /  is  the  number  of  rigid  body  nodes.  A  necesszuy  condition 
for  the  resulting  stiffness  matrix  to  have  sufficient  rank  is  m>n-/,  where  m  is  the  number 
of  independent  ^>stress  parameters.  Because  the  derivation  of  the  element  stiffness  matrix 
involves  inversion  of  [H]  matrix  which  is  of  the  order  mxm,  it  is  advantageous  to  use  as 
minimum  a  number  of  parameters  as  possible.  There  are  also  indications  that  overuse  of 
^-parameters  will  yield  over-rigid  elements.  An  ideal  situation  is  to  choose  the  stress  terms 
such  that  m  is  equal  to  n-/.  However,  it  is  not  always  possible  to  adhere  to  this  as  the 
number  of  fi-terms  is  also  governed  by  the  fact  that  the  assumed  stress  polynomial  should 
have  as  many  terms  as  required  to  satisfy  the  equilibrium  equation  (12).  Hence  for  some 
elements,  number  of  ^-terms  may  be  more  than  n-/.  The  resulting  element  properties  by  this 
approach  will  in  general  not  be  invariant.  However,  if  an  optimal  local  reference  coordinate 
system  is  used  in  the  element  formulation,  invariance  may  be  achieved. 

To  illustrate  the  assumed  stress  hybrid  method,  components  of  the  element  stiffness 
matrix  are  derived  as  follows. 

The  stress  functions,  involving  5  /3-terms,  that  satisfy  the  equation  of  equilibrium  may 


'  1  y  0  0  0 
0  0  1x0 
.00001 

This  is  suitable  for  elements  having  up  to  4  nodes.  However  for  elements  with  higher 
number  of  nodes,  as  in  the  case  of  voronoi  elements,  it  is  necessary  to  increase  the  number 
of  /?- terms  to  get  rank  sufficiency  of  the  stiffness  matrix.  The  stress  polynomial  [P]  for  10  ^ 
terms  is  given  as 


■  1 

.V 

0 

0 

0 

X 

0 

0 

x^ 

0 

0 

1 

X 

0 

0 

y 

0  x^ 

0 

0 

0 

0 

1 

-jr 

—X 
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be  assumed  as  follows 
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and  for  17  terms  as 

'  1  y  0  0  0  X  0  y^  0  xy  0  Zx^y  ty^  0 

001x0  OyO®^  y^  0  xy  Zxy^  0  y^  x^  x^y 

00001  ~y— xO  0  -2xy  -y^/2  -xV2  -3x*y  0  -3xy^  “-xV3 

The  stress  polynomial  involving  17  ;?-tenns  is  thus  smtable  for  elements  having  up  to  10 
nodes  (number  of  degrees  of  freedom  2x10,  minus  number  of  rigid  body  nodes  3,  dictates  a 
minimum  of  17  j5-tenns). Hence  for  elements  having  up  to  4  nodes,  5  ^-terms  are  used;  for 
elements  having  more  than  4  but  less  than  or  equal  to  6  nodes,  10  ^-terms  are  used;  for 
elements  having  more  than  6  but  less  than  or  equal  to  10  nodes,  17  terms  are  used. 

For  the  ith  side  of  an  element 


{«}  =  [LM 


(  ^2t-l  1 

1-a//,-  0  a/li  0  I  ^2t  I 

0  1  -  a/li  0  a/li  J  I  92<+i  | 

V  ?2t+2  J 


where  U  is  the  length  of  side  i  and  a  is  the  distance  measured  from  node  i  to  node  :+l.  The 
assumed  boundary  displacement  variation  is  linear.  Therefore,  the  interelement  compatibil¬ 
ity  is  assured  if  the  nodal  displacements  of  adjacent  elements  2ire  matched. 

Traction  on  the  side  i  is  given  as  (for  the  sake  of  simplicity  only  5  terms  are  considered) 


{T}i  = 


+  Tgyni2 
[  Txynn  +  a'yni2 


nil 

0 


t/nn 

0 


0  0  n,-2 

n,*2  xn,-2  n,i 


[Rm 


from  which 


nii(l-a/li) 

0  nua/U 

0 

yn,i(l  -  a//.) 

0  ynna/li 

0 

0 

'isCl  -  a/li)  0 

ni2a/li 

0 

xn.j(l  -  a/li)  0 

X7T|2  n  j  li 

n,j(l  -  a/li) 

nail -a/li)  n^a/U 

fill  a//,- 

and  [G]i  = 

The  contribution  to  components  G/,j 

for  the  side  i  is  given  by 

J/7  2:  -  1 

2i 

2i  +  1 

2z  +  2 

0 

n.xli/2 
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2  n,i/j(y,/3  +  j/,+i/6) 
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nuliivi/S  +  Vi+x/Z) 

0 
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nolij^ 
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It  is  interesting  to  note  that  for  a  four  node  element  if  only  three  ^-terms  are  used,  the 
element  sti&ess  matrix  obtained  by  assumed  stress  hybrid  method  becomes  exactly  same  as 
that  obtained  by  Mallett’s  finite  difference  element  scheme.  Usage  of  three  )9-terms  means 
assumption  of  constant  ffy  and  within  the  element  which  is  also  the  case  in  finite 
difference  element  scheme. 


3,2  Invariance  of  Stiflfhess  Matrix 

The  stifiEness  matrix  of  an  element  calculated  on  the  basis  of  global  coordinates  is  not  invari¬ 
ant  as  not  all  the  assumed  stresses  are  complete  polynomiads  of  the  same  degree.  However, 
invariance  may  be  achieved,  as  suggested  by  Cook  [3],  without  changing  the  stress  polynomial 
by  calculating  the  element  stiffness  in  a  local  coordinate  system  having  a  fixed  orientation 
with  respect  to  the  element,  regardless  of  how  the  element  may  be  oriented  in  global  coordi¬ 
nates.  Thus  if  a  given  structure  and  its  loads  are  rotated  some  amount  in  global  coordinates, 
nodal  loads  associated  with  each  element  maintain  their  original  direction  with  respect  to 
the  local  element  coordinates.  Thus  the  element  response  becomes  invariant  as  it  is  not 
affected  by  the  rotation. 

Once  a  local  coordinate  system  is  defined  the  calculation  proceeds  in  the  following  man¬ 
ner. 

1.  Nodal  coordinates  of  an  element  is  transferred  to  the  local  system  and  the  stiffness 
matrix  is  calculated  in  the  local  system. 

2.  Element  stiffness  matrix  in  the  local  system  is  then  transformed  to  the  global  coordinate 
system  and  assembled  in  the  global  stifness  matrix.  The  transformation  is  accomplished 
by  the  operation 

[k]global  =  [An[^]/oca/[A]  (23) 

where  [A]  is  the  coordinate  tranformation  matrix. 

3.  .A.fter  assembly,  element  nodal  displacements  are  calculated  in  the  global  xy  directions, 

4.  To  compute  stresses,  stress  parameters  ^  calculated  in  the  local  system  are  used. 
{(^}ioeai  =  [P]iocai{^} local  gi^cs  stesses  in  the  local  system  which  are  then  transformed 
to  the  global  system. 

Choice  of  finding  the  angle,  the  local  system  makes  with  global  system,  is  not  unique  and 
alternative  definitions  of  this  angle  are  possible.  However,  finding  a  local  coordinate  system 
for  a  polygonal  element  is  a  problem  as  a  polygon  does  not  exhibit  any  preferred  direction. 
In  the  present  formulation,  the  local  coordinate  of  an  element  is  defined  by  two  different 
methods. 
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In  the  first  method,  one  of  the  axes  of  the  local  coordinate  system  is  aligned  with  the 
longest  edge  of  the  element.  In  other  words,  the  angle  the  longest  edge  makes  with  one  of  the 
global  axes,  is  the  angle  by  which  the  global  axes  are  rotated  to  obtain  the  local  coordinate 
axes  for  the  particular  element. 

In  the  second  method,  an  element  is  first  mapped  to  a  master  4,  8  or  12  node  plane 
isoparametric  element  in  s-t  space.  If  the  element  has  less  than  or  equal  to  4  nodes,  it  is 
mapped  to  a  4  node  master  element,  if  it  has  nodes  between  5  and  8  it  is  mapped  to  an  8 
node  master  element  and  so  on.  Wherever  necessary,  extra  nodes  are  created  in  the  element 
mid  sides  so  that  there  is  a  one  to  one  correspondence  between  the  nodes  of  the  actual  and 
master  element.  Jacobian  matrix  for  this  transformation  is  then  calculated  and  decomposed 
into 


[J]  =  [R][U] 

« 

where  [R]  represents  a  pure  rotational  transformation  matrix  and 


(C  +  Vd^I) 


(24) 


where  [C]  = 

From  the  above  relationship  [i2]  is  obtained  as 

(fl)  =  [JlltT'  (25) 

This  rotational  transformation  matrix  then  is  used  for  rotating  the  global  axie  to  obtain 
the  local  coordinate  axes.  Rest  of  the  calculation  proceeds  as  outlined  earlier. 


3.3  Numerical  Examples 

Stretching  Problem: 

A  square  material  body  of  size  40  mm  x  40  mm  is  discretized  into  24  voronoi  polygons, 
as  shown  earlier  in  Fig.  1,  and  is  subjected  to  a  uniform  tensile  force  of  100  KN.  Analysis 
of  this  problem  by  the  Assumed  Stress  Hybrid  Method  yields  exact  end  deflection  of  4.55 
mm  in  X  direction  and  o’,  of  2.5  kN/sq.mm  in  ail  elements.  Also  the  displacement  of  the 
nodes  in  X  direction  versus  their  position  from  the  fixed  end  has  been  found  to  be  linearly 
proportional  (Fig.  2).  Thus  both  Mallett’s  method  and  assumed  stress  hybrid  method  yield 
exact  results  of  a  stretching  problem. 

The  material  body  is  now  aligned  at  45®  to  the  global  axes  and  is  subjected  to  the  same 
stretching  load  (Fig.  6).  The  results  of  deflection  and  stress  in  the  direction  the  load  are 
again  found  to  be  4.55  mm  and  2.5  KN/sq.mm  which  are  exact  figures. 
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Fig.  6.  Body  oriented  at  45®  and  subjected  to  tensile  load. 


Simple  Shear  Problem: 

The  same  40  mm  x  40  mm  material  body  discretised  with  the  same  24  voronoi  polygons 
is  subjected  to  suitable  displacement  boundary  condition  so  as  to  simulate  simple  shear  (Fig. 
7).  Assumed  stress  hybrid  method  yielded  exact  value  of  =  0.0385  KN/sq.mm  and  cr^  = 
cTy  =  0  in  every  element. 


Fig.  7.  Body  subjected  to  simple  shear. 
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Distiocc  iloof  y  dtrcctka 


Fig.  8.  Shear  stress  across  section  AA 


Bending  Problems: 

Example  1. 

A  material  body  of  size  40  mm  x  160  mm  is  used  as  a  cantilever  beam  with  one  end 
fixed  and  other  end  subjected  to  pure  moment  of  40  KN-mm.  The  beam  is  discretized 
into  28  voronoi  polygon  elements  (Fig.  9).  Longitudinal  stress  along  a  transverse  section 
calculated  by  assum^  stress  hybrid  method  conformed  more  closely  to  the  theoretical  values 
than  that  obtained  by  Mallett’s  method  (Fig.  10).  The  vertical  displacements  of  the  nodes 
at  the  top  surface  agree  very  well  with  the  theoretical  values  (Fig.  11).  The  calcui  ted  tip 
deflection  of  4.175  mm  has  an  error  of  4.4%  with  respect  to  the  theoretical  deflection  of  4.368 


mm. 


Slgna  X  al  ss3.7  cm 


Fig.  9.  Beam  subjected  to  pure  end  moment  (28  elements). 


* 


Fig.  10.  Sigma  X  along  a  transverse  section. 
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ak»f  X  dincliaii 


Fig.  11,  Deflection  of  top  surface. 


Examples  2  and  3. 


The  same  40  mm  x  160  mm  beam  is  now  discretized  into  53  and  152  voronoi  polygon 
elements  respectively  (Fig.  12  &  13)  and  subjected  to  piire  end  moment  of  40  KN-mm. 
<7,  along  a  transverse  section  is  plotted  (Fig.  14)  for  all  three  cases,  i.c.  beam  discretized 
by  28  elements,  53  elements  and  152  elements.  The  calculated  tip  deflection  for  the  beam 


with  53  elements  is  4.346  mm,  an  error  of  0.5%  compared  to  the  theoretioatdeflection.  The 
corresponding  figures  for  the  beam  with  152  elements  are  4.452  mm  an^^^^  respectively. 


Fig-  12.  Beam  subjected  to  pure  end  moment  (53  elements). 


Fig.  13.  Beam  subjected  to  pure  end  moment  (152  elements). 


Distance  along  Y  direction 


Fig.  14.  Sigma  X  along  a  transverse  section. 


Fig.  14  indicates  that  as  more  and  more  elements  are  used  to  discretize  the  beam,  the 
stress  distribution  converges  towards  the  theoretical  value.  Fig.  15  shows  the  distribution 
of  <7r  throughout  the  beam  discretized  by  152  elements.  Except  near  the  area  of  application 
of  moment  load,  where  the  Saint  Venant  effect  is  predominant,  the  (Tx  stress  distribution 
conforms  well  with  the  theoretical  stress  distribution. 
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Fig.  15.  Distribution  of  Sigma  X  in  beam  having  152  elements. 


The  beam,  discretized  into  53  elements,  is  now  aligned  at  45®  with  the  global  axis(Fig. 
16).  Stresses  in  some  elements  deteriorate  if  the  calculations  are  based  on  global  coordi¬ 
nate.  However,  use  of  a  local  coordinate  system  for  every  element  improves  the  results  and 
makes  the  elements  invariant  with  respect  to  the  orientation  of  the  beam  with  the  global 
axis  (Fig.  17).  Local  coordinate  for  an  element  is  usually  obtained  by  rotating  the  global 
system  until  one  of  the  axes  aligned  with  the  largest  edge  of  the  element,  and  is  identified  as 
Cook’s  method.  In  this  example  the  matrix  used  to  transform  the  global  system  to  the  local 
coordinate  system  is  also  obtained  by  decomposing  the  Jacobian  matrix,  as  detailed  earlier, 
and  is  identified  as  Jacobiein  method.  Results  based  on  local  coordinate  system,  obtained  by 
either  method,  compare  well  with  each  other. 
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l.oiif;iliulinal  Stress  at  distance  7.6  from  fixed  end 


X 

Fig.  16.  Beam  oriented  at  45®  and  subjected  to  end  moment(53  elements). 


Distance  aionc  transrerse  direction 

Fig.  17.  Longitudinal  stress  along  a  transverse  section  using  Local  &  Global  coprd. 


XHc  above  exercise  is  cdxried  out  on  the  beam  with  152  elements,  aligned  at  45  with  the 
global  axis.  Results  based  on  local  coordinate  system  are  found  to  improve  (Fig.  18).  ^  ^ 


Fig.  18.  Longitudinal  stress  along  a  transverse  section  using  Local  k  Global  coordinate. 

(152  elements) 


Example  4. 

Having  shown  that  the  results  improve  as  the  number  of  elements  increase,  the  same  40 
mm  X  160  mm  beam  is  discretized  into  86  elements  (Fig.  19),  most  of  which  are  of  regular 
hexagonal  shape,  and  the  same  exercise  is  earned  out.  Fig.  20  indicates  that  longitudinal 
stresses  in  the  elements  across  a  transverse  section  agree  very  well  with  the  theoretical  values. 
Fig.21  shows  the  distribution  of  Cg  in  the  entire  beam.  Except  near  the  area  of  application 
of  moment  load,  where  the  Saint  Venant  effect  is  predominant,  the  cTx  stress  distribution 
conforms  very  well  with  the  theoretical  stress  distribution. 

Results  of  this  example  indicate  that  as  the  elements  become  more  regular  in  shape,  the 
results  improve. 


21 


SicHW  X  at 


- -a-X 

Fig.  19.  Beam  subjected  to  pure  end  moment  (86  elements). 


Disunce  along  Y  direction 


Fig.  20.  Sigma  X  along  a  transverse  section  (86  elements). 
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Fig.  21.  Distribution  of  Sigma  X  in  beam  having  86  elements. 


The  beam  is  now  oreinted  at  45°  with  the  global  axis  (Fig.  22),  Fig.  23  indicates  that 
local  coordinate  system  for  every  element  produces  better  results  than  that  obtained  by  a 


Fig.  22.  Beam  oriented  at  45°  and  subjected  to  end  moment. 


Distance  along  trans>’erse  direction 


Fig.  23.  Longitudinal  stress  along  a  transverse  section  using  Local  &  Global  coord. 


Example  5. 

The  material  body  of  size  40  mm  x  160  mm,  with  regular  86  elements,  is  used  as  a 
cantilever  beam  subjected  to  a  downward  acting  load  of  0.5  KN  at  the  end  (Fig.  24).  The 
load  is  distributed  over  the  edge  in  a  parabolic  distribution  such  that  the  total  area  under 
the  parabola  is  equal  to  the  load  0.5  KN.  Shear  stress  across  a  transverse  section  is  plotted 
in  Fig.  25  and  found  to  agree  well  with  the  theoretical  distribution. 


Fig.  24.  Beam  under  end  load 
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Fig.  25.  Shear  stress  distribution  along  a  transverse  section. 

* 


25 


4  Modelling  of  Second  Phase  Material  Present  in  Voronoi 
Cell 

Having  established  a  finite  element  formulation  for  polygonal  elements  with  isotropic  material 
properties,  next  step  is  to  account  for  the  second  phase  particulate  present  in  every  element. 

The  problem  of  modelling  the  presence  of  a  second  phase  material  within  a  finite  element  was 
addressed  by  Accorsi  [4]  by  introducing  a  transformation  strain  in  those  regions.  The  finite 
element  method  reflects  the  shape,  size  and  location  of  the  second  phase  material  within  the 
element. 


4.1  Problem  Formulation 

The  formulation  is  based  on  dividing  the  total  problem  into  two  separate  subproblems. 
First,  called  homogeneous  problem,  corresponds  to  the  finite  element  analysis  of  the  body 
without  material  discontinuity  and  the  second,  called  the  deviation  problem  is  defined  by  the 
difference  between  the  actual  and  homogeneous  field  quantities.  The  solution  of  the  actual 
problem  can  be  found  by  adding  the  corresponding  solutions  of  homogenous  and  deviation 
problems. 

Let  superscripts  *<,  o  and  /  indicate  the  actual  or  nonhomogenous  problem,  the  homoge¬ 
nous  problem  and  the  deviation  problem  respectively. 

The  constitutive  relation  for  the  actual  material  is  thus  given  as 


(7*  =  C{x)€^  (26) 

where  C{x)  is  the  elasticity  matrix. 

The  transformation  strain  mentioned  earlier  is  defined  by  the  following  equation. 

(7*  =  axy  =  C-)  (27) 

where  is  the  elasticity  matrix  of  homogenous  material  and  e*  is  the  strain  in  the  discon¬ 
tinuity.  From  the  above  relationship  one  obtains 

c'  =  [C^  -  =  H{x)r€^  (2S) 

where  H{x)  =  [C^  -  C(x)]-' 

Now  as  per  earlier  definition 

or 

cV  e'  =  H(x)C\- 

The  constitutive  relation  for  the  delation  problem  is  given  as 

=  =  (31) 


(29) 

(30) 
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Discretization  of  the  homogenous  problem  is  straight  forward.  Using  the  principle  of 
minimum  potential  energy,  the  element  equations  are  obtained.  The  assembled  equations 
for  the  homogenous  problem  are  written  as 

=  F"  (32) 

where  6^,  F®  are  the  global  stiffness  matrix,  global  displacement  vector  and  global 
load  vector  respectively.  For  the  purpose  of  solving  the  above  equation,  they  can  be  rewritten 
as 


g.  =  (33) 

where  K^,Fm  indicate  iir®,F®  modified  for  displacement  boundary  conditions. 
Similarl}',  the  element  equation  for  the  discretized  deviation  problem  is  written  as 


Kx  =  f: 

where 

F:=  f  B^Cydv 

B  is  the  strain  displacement  matrix  such  that 


(34) 

(35) 


<  =  BS:  (36) 

and  F®  is  the  elemental  stiffness  matrix. 

For  an  element  of  matrix  material  C®  containing  me  microstructural  discontinuities  with 
material  property  =  l,2..me  and  occupying  regions  Hr,  the  transformation  strain  is 
non  zero  only  within  the  region  Hr-  The  equation  (35)  becomes 


me  - 

B^cydv 

rsri*'”’* 


(37) 


As  the  transformation  strain  £"  is  constant  within  a  discontinuity  the  above  equation  (37) 
is  further  reduced  to 


where 


b:  = 


Assembly  of  the  equation  (34)  gives 


K^S'  =  F- 


(38) 


(39) 
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where 


e 

<t>t  is  a  transformation  matrix  which  when  operated  on  global  displacement  vector  8 
returns  the  element  displacement  vector  8^. 

Equation  (39),  when  modified  for  displacement  boundary  condition,  takes  the  form 

K18'  =  rF-  (40) 

Since  the  displacement  boundary  conditions  for  the  homogenous  and  deviation  problems 
correspond  to  the  same  nodes,  the  modification  to  the  stiffness  matrix  is  same  for  both  the 
problems.  For  the  deviation  problem,  wherever  the  displacement  boundary  conditions  are 
zero,  zeros  are  placed  in  the  corresponding  positions  in  F*.  To  carry  out  this  modification 
of  F*  a  matrix  ^  is  used.  is  an  identity  matrix  with  zeros  in  the  diagonal  in  the  positions 
corresponding  to  the  displacement  boxmdary  conditions. 

From  equation  (40) 

F=(F»)-VF-  (41) 

Therefore 

K  =  OtF 


_  me 

e 

=  o,{Klr^rbY  (42) 

SSZl 

where  M  is  the  total  number  of  discontinuities  in  the  system. 

Writing  equation  (30)  in  terms  of  nodal  displacements  for  an  element,  one  obtains 

BS:  +  =  H{x)cy  (43) 

Integrating  equation  (43)  over  fir,  the  domain  of  rth  discontinuity  in  the  element,  one 
obtains 

RlSl  +  =  VrH^cy^ 
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where  K  is  volume  of  rth  discontinuity  and  =  [C®  — 

Rearranging  the  above  equation 

VrH:Cy^  -  RX  =  Ris:  =  (44) 

Substituting  the  value  of  (5'  from  equation  (42)  into  equation  (44)  we  have 

£  [v.H:S^‘  -  cy’  =  Rloy  (45) 

where  6^*  is  the  Kronecker  delta.  The  above  equation  is  now  assembled  by  evaluating  the 
terms  for  each  discontinuity  in  the  system  r.s  =  The  assembled  equation  takes 

the  form 

5(C"£")  =  R6^  (46) 

It  may  be  observed  from  equation  (45  )  that  e*’’  is  function  of  i?;  and  ifj.  i^e(=  Bdv) 
in  turn  is  dependent  upon  the  shape  and  size  of  the  discontinuity,  as  the  integration  is  over 
Qr?  and  also  on  the  location  of  the  discontinuity  within  the  element  due  to  the  presence  of 
B  matrix  in  the  expression.  If  the  element  to  start  with  is  not  constant  strain  type,  B  is  a 
function  of  the  coordinates  within  the  element.  matrix  contains  the  material  property  of 
the  discontinuity.  Thus  c*’*  depends  upon  the  shape,  size  and  location  of  the  discontinuity 
within  the  element. 

The  subroutine  for  calculation  of  6^  follows  the  following  steps 
-From  equation  (46)  cdculate 

-Calculate  by  substituting  the  value  of 

* 

-Find  global  F-  =  Ee 

-Substitute  F*  in  equation  (41),  to  get  F  = 

A  standard  2-D  finite  element  program  FEM-2  [5]  is  used  to  calculate  and  i.e. 
the  solution  of  the  homogenous  and  deviation  problem  respectively.  Solution  of  the  actual 
problem  is  thus  +  F. 
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4.2  Numerical  Examples 

Stretching  Problem: 

The  homogenization  formulation  allows  the  discontinuity  to  be  treated  as  a  void  by  set¬ 
ting  its  material  properties  to  zero.  A  square  material  body  with  a  central  hole  is  subjected 
to  tensile  load.  Resulting  stress  distribution  along  a  transverse  section  passing  through  the 
hole  is  calculated  by  conventional  FEM  and  FEM  coupled  with  homogenization  method.  For 
the  first  case  (Fig.  26)  material  aroimd  the  hole  is  discretized  by  four  node  elements  and  for 
the  second  case  (Fig.  27),  the  material  is  discretized  in  such  a  way  that  the  entire  hole  is 
contained  within  an  element. 
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Stress  distribution  (Fig.  28)  by  either  method  compare  well  up  to  a  distance  j  from  the 
edge  of  the  hole,  D  being  the  hole  diameter.  After  that,  the  sharp  rise  in  stress  as  the  edge  is 
approached  is  not  reflected  properly  by  homogenization  method.  This  however,  is  expected 
as  the  homogenization  method  models  the  presence  of  the  entire  hole  in  an  element  in  an 
average  sense. 


Fig.  28.  Sigma  X  across  a  section  through  the  hole. 


Bending  Problem: 

A  cantilever  beam,  having  three  second  phase  materials  and  subjected  to  a  pure  mo¬ 
ment  at  the  end  is  analyzed  by  conventional  FEM  and  FEM  coupled  with  homogenization 
method.  For  the  mventional  FEM  the  beam  is  discretized  as  shown  in  Fig.  29  and  for  the 
homogenization  method  the  discretization  is  shown  in  Fig.  30.  The  Young^s  modulus  of  the 
second  phase  material  is  assumed  as  4  KN/sq.mm  which  is  double  that  of  matrix  material. 
Poisson’s  ratio  is  assximed  as  0.3. 
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Fig.  29.  Analysed  by  FEM2 


- ►! 

Fig.  30.  Analysed  by  FEM2  with  Homogenization. 


The  stress  distribution,  along  a  transverse  section  passing  through  the  second  phase 
materials,  obtained  by  either  method  is  plotted  in  Fig.  31.  The  stress  predicted  by  homoge¬ 
nization  techniciue  is  found  to  compare  somewhat  well  in  the  matrix  material,  however  wi  thin 
the  second  phase  material,  the  stress  is  well  above  than  that  predicted  by  conventional  FEM. 
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Fig.  31.  Sigma  X  across  a  section  through  the  second  phase  material. 

* 

All  the  numerical  results  of  conventional  FEM  and  homogenization  methods,  discussed  so 
far,  are  based  on  a  program  for  analysis  of  linearly  elastic  structures  in  a  plane  using  4-node 
isoparametric  elements,  called  FEM2.  However,  for  the  analysis  of  composites  with  randomlj 
dispersed  second  phase  materials,  the  elements  that  will  be  encountered,  as  explained  earlier 
in  the  report,  will  have  more  than  four  nodes.  For  the  analysis  of  these  polygonal  elements, 
assumed  stress  hybrid  method  is  introduced.  The  present  problem  of  a  beam  dispersed  with 
three  second  phase  materials  and  subjected  to  a  pure  end  moment  is  then  analysed  by  the 
homogeneation  technique  coupled  with  assumed  stress  hybrid  method.  The  results,  as  shown 
in  Fig.  30  compare  very  well  with  that  obtained  by  homogenization  method  coupled  with 
7EM2. 

The  same  beam  bending  problem  is  solved  by  discretizing  the  beam  into  irregular  polygo¬ 
nal  elements  (Fig.  32)  and  using  homogenization  coupled  with  assumed  stress  hybrid  method. 
Longitudinal  stresses  along  the  same  transverse  section  using  irregular  and  regular  mesh  are 
plotted  in  Fig.  33.  The  results  are  found  to  compare  well 
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Fig.  32.  Beam  discretized  by  polygonal  elements. 
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Fig.  33.  Sigma  X  across  a  section  through  the  second  phase  material. 

Finally  the  beam  of  size  40  mmx  160  mm  which  was  used  in  several  earlier  examples 
analyst  by  homogenization  coupled  with  assumed  stress  hybrid  method.  The  beam 


Sigma  X  at  aslJJ 


however,  is  now  dispersed  with  64  second  phase  particulates  and  is  subjected  to  a  pure  end 
moment  of  40  KN-mm(Fi^.  34).  Longitudinal  stresses  along  a  transverse  section  which 
passes  through  several  second  phase  particulates  are  plotted  in  Fig.  35. 


t - ►x 

Fig.  34.  Composite  beam  subjected  to  pure  end  moment  (64  elements) 


Fig.  35.  Sigma  X  across  a  section 
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